#include <iostream>
#include "spline.h"
#include <cmath>
#include <vector>
#include <algorithm>
#include <fstream>
class Fun : public Function
{
public:
    double operator()(double _x)
    {
        return 1.0/(1.0 + 25.0*_x*_x);
    }
    double diff(double _x)
    {
        return -50.0*_x/pow(1 + 25.0*_x*_x, 2);
    }
    double diff_2(double _x)
    {
        return 50.0*(75.0*_x*_x - 1)/pow(1.0 + 25.0*_x*_x, 3);
    }
};

int main()
{
    std::ofstream fout_point, fout_error;
    fout_point.open("A_point.txt");
    fout_error.open("A_error.txt");
    int N[8] = {6, 11, 21, 41, 81, 161, 321, 641};
    int n = 8;
    double a = -1;
    double b = 1;
    int order = 3;
    int way[3] = {1, 2, 3};
    Fun f;
    for(int k = 0; k <= 2; k++)
    {
        for(int i = 0; i <= n-1; i++)
        {
            std::vector<double> Point;
            for(int j = 0; j <= N[i]-1; j++)
            {
                Point.push_back(a + j*(b-a)/(N[i]-1));
            }
            ppForm_interpolation P(f, Point, way[k], order);//ppform
            P.solve();
            for(int m = 0; m < 100; m++)
            {
                fout_point << P(-1 + m*2.0/99) << std::endl;
            }
            double max_error = 0;
            for(int j = 0; j < N[i]-1; j++)
            {
                double error = fabs(P(0.5*(Point[j]+Point[j+1])) - f(0.5*(Point[j]+Point[j+1])));
                if(error > max_error)
                {
                    max_error = error; 
                }
            }
            fout_error << max_error << std::endl;
        }
        for(int i = 0; i <= n-1; i++)
        {
            std::vector<double> Point;
            for(int j = 0; j <= N[i]-1; j++)
            {
                Point.push_back(a + j*(b-a)/(N[i]-1));
            }
            B_spline_interpolation B(f, Point, way[k], order);//B_spline
            B.solve();
            for(int m = 0; m < 100; m++)
            {
                fout_point << B(-1 + m*0.02) << std::endl;
            }
            double max_error = 0;
            for(int j = 0; j < N[i]-1; j++)
            {
                double error = fabs(B(0.5*(Point[j]+Point[j+1])) - f(0.5*(Point[j]+Point[j+1])));
                if(error > max_error)
                {
                    max_error = error; 
                }
            }
            fout_error << max_error << std::endl;
        }
    }
    fout_point.close();
    fout_error.close();
    return 0;
}
